//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ephi.hpp"

class myElement : public StaticElement
{
public:
  myElement (const vect3d& pos, prec_t R, prec_t Q) : pos(pos), R(R), Q(Q) { }

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield)
  {
    vect3d r = pos - this->pos;
    prec_t r2 = r*r;
    prec_t rlen = prec_t_sqrt(r2);
    prec_t r3 = r2 * rlen;

    if ( rlen < 1e-10 )
      {
        return;
      }
    else if ( rlen >= R )
      {
        efield += (Q / r3) * r;
      }
    else
      {
        prec_t r5 = r3 * r2;
        prec_t r7 = r5 * r2;
        prec_t R2 = R * R;
        prec_t R4 = R2 * R2;
        prec_t R5 = R4 * R;
        prec_t R7 = R5 * R2;
        prec_t R9 = R7 * R2;
        efield += (Q * (35 * r3 / (8 * R5) - 21 * r5 / (4 * R7) + 15 * r7 / (8*R9))) * r.normal();
      }
  }

  void addPotential (const vect3d& pos, prec_t& potential)
  {
    vect3d r = pos - this->pos;
    prec_t r2 = r*r;
    prec_t rlen = prec_t_sqrt(r2);
    prec_t r4 = r2 * r2;
    prec_t r6 = r4 * r2;
    prec_t r8 = r6 * r2;
    if ( rlen >= R )
      {
        potential += Q / rlen;
      }
    else
      {
        prec_t R2 = R * R;
        prec_t R4 = R2 * R2;
        prec_t R5 = R4 * R;
        prec_t R7 = R5 * R2;
        prec_t R9 = R7 * R2;

        potential += -Q * (35 * r4 / (32*R5) - 21 * r6 / (24 * R7) + 15 * r8 / (64 * R9) - 93 / (64 * R));
      }
  }

  void draw(Screen& screen, const Screen::color& c)
  {
  }

  std::string get_povray_script ()
  {
    return std::string();
  }

  bool inContact (const vect3d& pos)
  {
    return false;
  }

private:
  const vect3d& pos;
  prec_t R;
  prec_t Q;
};


int main (int argc, char *argv[])
{
  Statics statics;
  statics.addStaticElement (new myElement(vect3d(0.0, 0, 0), 0.5, 1));
  ElectroDynamics dyn(statics);
  Scene scene(dyn);
  Screen screen(800, 800, 2, 2);
  vect3d bf, ef;
  statics.getFields (vect3d(0.0, 0, 0), bf, ef);
  bf.print("bf");
  ef.print("ef");
  scene.calc (screen, Scene::CALC_EFIELD);
  scene.set_coloring (0, false);
  scene.render_map (screen);
  screen.write ("tt.bmp");
  scene.calc (screen, Scene::CALC_POTENTIAL);
  scene.set_coloring (32, false);
  scene.render_map (screen);
  screen.write ("tt_pot.bmp");
#if 0
  PTriangle t(vect3d(1, 0, 0), vect3d(3, 0, 0), vect3d(2, 1, 1));
  prec_t pot = 0;
  t.addPotential (vect3d(2, 0.5, 2), pot);
#endif
#if 0
  Screen screen(800, 800, 2.0, 2.0);
  Statics statics;
  ElectroDynamics dyn(statics), dyn2(statics);
  dyn.setDTM(100, 1);

  statics.addStaticElement (new StaticRing(5000, vect3d(0, 0, 0), 0.25, 0, vect3d(0.0, 0.0, 1.0)));
  statics.draw (screen);

  prec_t speed = ElectroDynamics::ms2eV (10000000, ELECTRON_MASS);
  dyn.inject_electron (vect3d(-0.5, 0, 0), vect3d(1, 0, 0), speed);

  vect3d pos = dyn.get_pos(0);
  vect3d v = dyn.get_v(0);
  vect3d bf;
  dyn2.getBField (pos, bf);

  prec_t q = -ELEMENTARY_CHARGE;
  prec_t m = ELECTRON_MASS;
  vect3d f = -q * (v % bf);
  prec_t rlen = m * v * v / f.length();
  vect3d r = -v % bf; // - for the electron charge
  r = r.normal() * rlen;
  vect3d center = pos + r;

  for ( size_t i = 0; i < 100000; i++ )
    {
      dyn.step();
      //dyn.get_pos(0).print();
      screen.set_pixel (dyn.get_pos(0), Screen::RED);
    }

  r.print("r");

  prec_t phi = 359.9 * 2.0 * PREC_PI / 360.0;
  prec_t cphi = cos(phi);
  prec_t sphi = sin(phi);
  for ( size_t i = 0; i < 100000; i++ )
    {
      vect3d normal(0, 0, 1);
      vect3d next = center - cphi * r - sphi * (normal % r);
      pos = next;
      r = center - pos;

      vect3d bf;
      dyn2.getBField (pos, bf);
      vect3d f = -q * (v % bf);
      prec_t rlen = m * v * v / f.length();
      r = r.normal() * rlen;

      center = pos + r;
      screen.set_pixel (pos, Screen::GREEN);
      //screen.draw_cross (pos, Screen::GREEN, 0.02);
      //screen.draw_cross (pos + r, Screen::GREEN, 0.02);
    }

  screen.draw_axes();

  screen.write ("tt.bmp");

  return 0;
#endif
}

